source('helpers.R')
library(MaxPac)
library(effectsize)
theme_set(theme_minimal())
suppressWarnings(dir.create('PaperFigures'))
We read data from the RetroDuration questionnaire. Translate responses using QTranslate() (google sheet). Read response to the initial_RetrospectiveDuration question. Extract estimated duration since last login as ##:##:## –> Hours:Minutes:Seconds. Compute ClockDuration as the difference between current time and last login.
Compute :
load(file.path(dirBlursday,'Daily_Login_Times.RData'))
# Extracting data:
#
# RetroDuration <- gimmeRdata(dirBlursday, UniqueName = 'RetroDuration') %>%
# filter(!Country %in% c('UK','CO', 'US')) %>%
# QTranslate() %>%
# filter(Question_Key %in% c('initial_RetrospectiveDuration')) %>%
# mutate(Date = lubridate::date(Local_Date)) %>%
# pivot_wider(id_cols = c('PID', 'Session','Run', 'Country', 'Local_Date'), names_from = Question_Key, values_from = Response) %>%
# # take only the first response if there are several ones
# group_by(PID,Session,Run,Country) %>%
# slice(1) %>%
# add_Demographics() %>%
# # For EstimatedDuration
# # Using 3 AM to change day
# mutate(Date = lubridate::as_date(Local_Date - lubridate::duration(hours = 3))) %>%
# left_join(select(daily_login_times, PID, Session, Local_Date, Date), by = c('PID','Date'), suffix = c('','_First_Login')) %>%
# mutate(Date = lubridate::as_date(Local_Date)) %>%
# mutate(ClockDuration = Local_Date - Local_Date_First_Login) %>%
# select(Country, PID, Session, Run, Local_Date, ClockDuration, initial_RetrospectiveDuration,
# Handedness, Sex, Age, Date) %>%
# mutate(M = str_match(initial_RetrospectiveDuration,
# '^.*?(\\d+).*?[::].*?(\\d+).*?[::].*?(\\d+).*?$'),
# h = suppressWarnings(as.numeric(M[,2])),
# m = suppressWarnings(as.numeric(M[,3])),
# s = suppressWarnings(as.numeric(M[,4])),
# EstimatedDuration = lubridate::duration(hour = h, minute = m, second = s),
# ClockDuration = lubridate::as.duration(ClockDuration),
# rEstimatedDuration = EstimatedDuration / ClockDuration,
# EstimationError = EstimatedDuration - ClockDuration,
# rEstimationError = EstimationError / ClockDuration) %>%
# select(-M) %>% # ,-h,-m,-s
# ungroup()
#
# save(RetroDuration, file = file.path(dirBlursday,'RetroDuration.RData'))
load(file = file.path(dirBlursday,'RetroDuration.RData'))
RetroDuration %>%
mutate(across(c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError), as.numeric)) %>%
pivot_longer(cols = c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError)) %>%
ggplot(aes(x = value, fill = Country)) +
geom_histogram(bins = 50) +
facet_wrap(~name ,scales = 'free') +
ggtitle('Measure distributions before outliers detection')
## Warning: Removed 1464 rows containing non-finite values (stat_bin).
Distributions of ClockDuration, EstimatedDuration, rEstimationError, and rEstimatedDuration are highly skewed to the left. First capping EstimatedDuration and ClockDuration to a decent duration, then removing 5% most extreme rEstimationError.
Rationale:
ClockDuration cannot be below 60 s: no questionnaire/task can be completed under this duration, and we want at least one.ClockDuration cannot be above 18000 s = 5 h: subjects cannot be focused for 5h in a row.EstimatedDuration, but 5 fold more permissive (we allow EstimatedDuration up to 5 * 18000s, and below 60 / 5 s).rEstimationError is systematically discarded. This is done for each Session and Country separately.We count these
# capped_outliers:
Duration_max <- 18000
Duration_min <- 60
Ratio_rejection_ClockDuration_to_EstimatedDuration <- 5
# trim_outliers:
probs <- c(.025, .975)
# RetroDuration_clean <- RetroDuration %>%
# filter(!is.na(rEstimatedDuration), rEstimatedDuration != 0)
# nrowsbefore_outliers <- nrow(RetroDuration_clean)
# RetroDuration_clean <- RetroDuration_clean %>%
# group_by(Country,Session) %>%
# mutate(outliers = ClockDuration > Duration_max | ClockDuration < Duration_min | EstimatedDuration > Duration_max * Ratio_rejection_ClockDuration_to_EstimatedDuration | EstimatedDuration < Duration_min / Ratio_rejection_ClockDuration_to_EstimatedDuration) %>%
# filter(!outliers) %>%
# mutate(outliers = find.outliers(rEstimationError,meth = 'trim', probs = probs)) %>%
# filter(!outliers) %>%
# ungroup() %>% select(-outliers)
# nrowsafter_outliers <- nrow(RetroDuration_clean)
#
# cat('Global proportion of outliers', prop_outliers_global <- 1 - nrowsafter_outliers / nrowsbefore_outliers)
#
# save(RetroDuration_clean, file = file.path(dirBlursday,'RetroDuration_clean.RData'))
load(file.path(dirBlursday,'RetroDuration_clean.RData'))
RetroDuration_clean %>%
group_by(Country,Session) %>%
mutate(across(c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError), as.numeric)) %>%
pivot_longer(cols = c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError)) %>%
ggplot(aes(x = value, fill = Country)) +
geom_histogram(bins = 50) +
facet_wrap(~name ,scales = 'free') +
ggtitle('Measure distributions after outliers detection')
RetroDuration_clean %>% group_by(Country, Session) %>%
summarize(n = n())
load(file = file.path(dirBlursday, 'RetroDuration_clean.RData'))
tostat <- RetroDuration_clean %>%
group_by(Country, PID, Session) %>%
select(Country,PID, Session, Local_Date,EstimatedDuration, rEstimatedDuration, EstimationError, Age, Sex, Handedness, ClockDuration) %>%
filter(Session %in% c('S1','SC')) %>%
summarize(Local_Date = unique(Local_Date),
EstimationError = mean(EstimationError) / 60,
ClockDuration = mean(ClockDuration) / 60,
EstimatedDuration = mean(EstimatedDuration) / 60,
rEstimatedDuration = mean(rEstimatedDuration),
rEstimationError = mean(EstimationError/ClockDuration),
Age = unique(Age)) %>%
mutate(logrEstimatedDuration = log10(rEstimatedDuration)) %>%
add_SubjectiveConfinementIndices() %>%
add_StringencyIndex() %>%
add_TimeOfDay()
N <- tostat %>%
group_by(Session) %>%
summarize(N = n())
m <- lm(log10(EstimatedDuration) ~ log10(ClockDuration) * Session, data = tostat)
plot(m)
summary(m)
##
## Call:
## lm(formula = log10(EstimatedDuration) ~ log10(ClockDuration) *
## Session, data = tostat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84487 -0.10233 0.01158 0.11975 1.10009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.12709 0.02891 4.396 1.17e-05 ***
## log10(ClockDuration) 0.89451 0.01938 46.159 < 2e-16 ***
## SessionSC 0.16688 0.06903 2.418 0.0157 *
## log10(ClockDuration):SessionSC -0.12640 0.04594 -2.751 0.0060 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1804 on 1739 degrees of freedom
## Multiple R-squared: 0.587, Adjusted R-squared: 0.5863
## F-statistic: 824 on 3 and 1739 DF, p-value: < 2.2e-16
anova(m)
effectsize::eta_squared(m)
confint(m)
## 2.5 % 97.5 %
## (Intercept) 0.07038706 0.18379810
## log10(ClockDuration) 0.85650572 0.93252236
## SessionSC 0.03149689 0.30227052
## log10(ClockDuration):SessionSC -0.21650737 -0.03629085
(fig1a <- emmeans(m,~ClockDuration+Session,at = list(ClockDuration = seq(min(tostat$ClockDuration,na.rm = F), max(tostat$ClockDuration,na.rm = F), length.out = 20))) %>%
summary(type = 'response') %>%
ggplot(aes(x = ClockDuration, y = response, ymin = lower.CL, ymax = upper.CL, col = Session, group = Session)) +
geom_abline(slope = 1) +
geom_point(aes(y= EstimatedDuration, ymin = NA, ymax = NA),alpha = .5, size = 1, data = tostat) +
geom_line(size = 1) +
geom_ribbon(col = NA, fill = 'grey',alpha = .2) +
scale_x_log10() +
scale_y_log10() +
ylab('Duration Estimate (min)') +
xlab('Clock Duration (min)') +
scale_fill_manual(values = Palette_SessionS1SC) +
scale_color_manual(values = Palette_SessionS1SC))
## Warning: Ignoring unknown aesthetics: ymin, ymax
summary(m)
##
## Call:
## lm(formula = log10(EstimatedDuration) ~ log10(ClockDuration) *
## Session, data = tostat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84487 -0.10233 0.01158 0.11975 1.10009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.12709 0.02891 4.396 1.17e-05 ***
## log10(ClockDuration) 0.89451 0.01938 46.159 < 2e-16 ***
## SessionSC 0.16688 0.06903 2.418 0.0157 *
## log10(ClockDuration):SessionSC -0.12640 0.04594 -2.751 0.0060 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1804 on 1739 degrees of freedom
## Multiple R-squared: 0.587, Adjusted R-squared: 0.5863
## F-statistic: 824 on 3 and 1739 DF, p-value: < 2.2e-16
summary(emtrends(m,~Session, 'log10(ClockDuration)', infer = T, offset = -1, adjust = 'tukey'))
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
anova(m)
qqnorm(residuals(m))
effectsize::eta_squared(m)
summary(emmeans(m,~ClockDuration+Session,at = list(ClockDuration = c(10,30,60))),type = 'response') %>%
mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
SE = SE * 60,
response = lubridate::as.duration(response * 60))
summary(emmeans(m,~ClockDuration+Session,at = list(ClockDuration = seq(15,20, by = .1))),type = 'response') %>%
mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
response = lubridate::as.duration(response * 60),
identity = response - ClockDuration) %>%
group_by(Session) %>%
filter(abs(identity) == min(abs(identity))) %>%
select(-response)
summary(emmeans(m,~ClockDuration, at = list(ClockDuration = seq(15,20, by = .1))),type = 'response') %>%
mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
response = lubridate::as.duration(response * 60),
identity = response - ClockDuration) %>%
filter(abs(identity) == min(abs(identity))) %>%
select(-response)
## NOTE: Results may be misleading due to involvement in interactions
tostat %>% group_by(Session) %>% summarize(mean_ClockDuration = mean(ClockDuration))
tostat %>%
ggplot(aes(x=Age,y = after_stat(ncount), fill = Session)) +
geom_histogram(position = 'dodge') +
scale_fill_manual(values = Palette_SessionS1SC) +
scale_color_manual(values = Palette_SessionS1SC) +
ylab('normalized count')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 45 rows containing non-finite values (stat_bin).
We first insert all covariates in the data, looking at all sessions, now.
load(file = file.path(dirBlursday, 'RetroDuration_clean.RData'))
tostat <- RetroDuration_clean %>%
group_by(Country, PID, Session) %>%
select(Country,PID, Session, Local_Date, EstimatedDuration, rEstimatedDuration, EstimationError, Age, Sex, Handedness, ClockDuration) %>%
filter(!is.na(rEstimatedDuration), rEstimatedDuration != 0) %>%
summarize(n = n(),
Local_Date = unique(Local_Date),
EstimatedDuration = mean(EstimatedDuration),
rEstimatedDuration = mean(rEstimatedDuration),
rEstimationError = mean(EstimationError/ClockDuration),
EstimationError = mean(EstimationError),
ClockDuration = mean(ClockDuration),
Age = unique(Age)) %>%
mutate(logrEstimatedDuration = log10(rEstimatedDuration)) %>%
add_SubjectiveConfinementDuration() %>%
add_SubjectiveConfinementIndices() %>%
add_Mobility() %>%
add_StringencyIndex() %>%
add_TimeOfDay()
## Joining, by = c("Country", "PID")
(N <- tostat %>%
group_by(Session) %>%
summarize(N = n()))
Subjects were repeatedly tested in S1-S4, and a different pool of subjects were tested in SC. We thus try to see if the error related to subjects in the intercept of EstimatedDuration is worth modeling as a random effect. The standard deviation of the estimated random effects is much smaller than that of the residuals, suggesting that random effects can be ignored in this case.
m <- lmerTest::lmer(log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement + Age + poly(Hour_Of_Day,3) + (1 | PID), data = tostat )
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement +
## Age + poly(Hour_Of_Day, 3) + (1 | PID)
## Data: tostat
##
## REML criterion at convergence: -843.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5107 -0.5250 0.0382 0.6064 8.3086
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 0.003692 0.06076
## Residual 0.035068 0.18726
## Number of obs: 2147, groups: PID, 1736
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.513e-02 2.882e-02 1.853e+03 1.219 0.22304
## Stringency_Index -6.230e-04 2.310e-04 1.699e+03 -2.697 0.00706 **
## Subjective_Confinement -7.860e-04 1.309e-03 1.935e+03 -0.600 0.54829
## Age 5.869e-05 2.802e-04 1.597e+03 0.209 0.83413
## poly(Hour_Of_Day, 3)1 -1.730e-01 2.089e-01 2.126e+03 -0.828 0.40760
## poly(Hour_Of_Day, 3)2 -3.149e-01 2.055e-01 2.139e+03 -1.532 0.12560
## poly(Hour_Of_Day, 3)3 3.172e-01 2.070e-01 2.139e+03 1.533 0.12553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Strn_I Sbjc_C Age p(H_O_D,3)1 p(H_O_D,3)2
## Strngncy_In -0.662
## Sbjctv_Cnfn -0.691 0.088
## Age -0.361 0.017 -0.027
## p(H_O_D,3)1 -0.048 0.014 0.024 0.059
## p(H_O_D,3)2 -0.061 0.065 0.028 0.004 0.011
## p(H_O_D,3)3 -0.003 0.051 -0.014 -0.044 -0.040 0.011
anova(m)
effectsize::eta_squared(m)
Indeed, running a simple lm on the data with the same fixed effects retrieves the same effects.
m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
##
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement +
## Age + poly(Hour_Of_Day, 3), data = tostat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2366 -0.1115 0.0067 0.1270 1.7101
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.374e-02 2.835e-02 1.543 0.12300
## Stringency_Index -6.727e-04 2.259e-04 -2.978 0.00293 **
## Subjective_Confinement -9.998e-04 1.292e-03 -0.774 0.43919
## Age 2.528e-05 2.734e-04 0.092 0.92634
## poly(Hour_Of_Day, 3)1 -1.439e-01 2.083e-01 -0.691 0.48972
## poly(Hour_Of_Day, 3)2 -2.736e-01 2.055e-01 -1.332 0.18307
## poly(Hour_Of_Day, 3)3 3.098e-01 2.070e-01 1.497 0.13463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.197 on 2140 degrees of freedom
## (267 observations deleted due to missingness)
## Multiple R-squared: 0.006255, Adjusted R-squared: 0.003469
## F-statistic: 2.245 on 6 and 2140 DF, p-value: 0.03654
anova(m)
effectsize::eta_squared(m)
plot(m)
Same reasoning with Mobility_Transit:
m <- lmerTest::lmer(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age + poly(Hour_Of_Day,3) + (1 | PID), data = tostat )
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit +
## Subjective_Confinement + Age + poly(Hour_Of_Day, 3) + (1 | PID)
## Data: tostat
##
## REML criterion at convergence: -836.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3938 -0.5133 0.0453 0.6010 8.4372
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 0.003564 0.0597
## Residual 0.035079 0.1873
## Number of obs: 2147, groups: PID, 1736
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.536e-02 2.972e-02 1.685e+03 1.863 0.062654 .
## Stringency_Index -1.408e-03 3.726e-04 1.417e+03 -3.780 0.000163 ***
## Mobility_Transit -7.562e-04 2.824e-04 2.091e+03 -2.678 0.007464 **
## Subjective_Confinement -7.120e-04 1.307e-03 1.932e+03 -0.545 0.585982
## Age -1.624e-05 2.810e-04 1.584e+03 -0.058 0.953914
## poly(Hour_Of_Day, 3)1 -1.815e-01 2.086e-01 2.125e+03 -0.870 0.384345
## poly(Hour_Of_Day, 3)2 -3.008e-01 2.052e-01 2.138e+03 -1.466 0.142868
## poly(Hour_Of_Day, 3)3 3.051e-01 2.067e-01 2.138e+03 1.476 0.140090
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Strn_I Mblt_T Sbjc_C Age p(H_O_D,3)1 p(H_O_D,3)2
## Strngncy_In -0.594
## Mblty_Trnst -0.251 0.786
## Sbjctv_Cnfn -0.663 0.037 -0.023
## Age -0.373 0.088 0.098 -0.029
## p(H_O_D,3)1 -0.050 0.022 0.017 0.023 0.060
## p(H_O_D,3)2 -0.053 0.022 -0.023 0.029 0.002 0.011
## p(H_O_D,3)3 -0.008 0.049 0.021 -0.015 -0.042 -0.039 0.011
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
anova(m)
effectsize::eta_squared(m)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
Indeed, running a simple lm on the data with the same fixed effects retrieves the same effects.
m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
##
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit +
## Subjective_Confinement + Age + poly(Hour_Of_Day, 3), data = tostat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21370 -0.10848 0.00998 0.12468 1.73176
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.240e-02 2.909e-02 2.145 0.03208 *
## Stringency_Index -1.457e-03 3.616e-04 -4.029 5.8e-05 ***
## Mobility_Transit -7.790e-04 2.808e-04 -2.774 0.00558 **
## Subjective_Confinement -9.164e-04 1.291e-03 -0.710 0.47772
## Age -4.676e-05 2.742e-04 -0.171 0.86461
## poly(Hour_Of_Day, 3)1 -1.559e-01 2.080e-01 -0.750 0.45354
## poly(Hour_Of_Day, 3)2 -2.607e-01 2.052e-01 -1.271 0.20402
## poly(Hour_Of_Day, 3)3 2.979e-01 2.067e-01 1.441 0.14965
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1967 on 2139 degrees of freedom
## (267 observations deleted due to missingness)
## Multiple R-squared: 0.009818, Adjusted R-squared: 0.006578
## F-statistic: 3.03 on 7 and 2139 DF, p-value: 0.003583
anova(m)
effectsize::eta_squared(m)
confint(m)
## 2.5 % 97.5 %
## (Intercept) 0.005345989 0.1194542821
## Stringency_Index -0.002166144 -0.0007477780
## Mobility_Transit -0.001329668 -0.0002283539
## Subjective_Confinement -0.003447168 0.0016144056
## Age -0.000584506 0.0004909846
## poly(Hour_Of_Day, 3)1 -0.563851130 0.2519801769
## poly(Hour_Of_Day, 3)2 -0.663098509 0.1416838641
## poly(Hour_Of_Day, 3)3 -0.107437203 0.7032299584
plot(m)
### Illustrating
m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
em <- summary(emmeans(m, ~ Stringency_Index, at = list(Stringency_Index = seq(20,100,by = 20))), type = 'response')
(fig1b <- ggplot(aes(x = Stringency_Index,y = rEstimatedDuration, col = Session), data = tostat) +
geom_jitter(alpha = .3, size = 1) +
geom_line(aes(y=response), col = 'black', data = em) +
geom_ribbon(aes(y=response, ymin = lower.CL, ymax = upper.CL), col = NA, alpha = .2, data = em) +
scale_y_log10( limits = c(.3,3)) +
scale_fill_manual(values = Palette_Session) +
scale_color_manual(values = Palette_Session) +
ylab('Relative Duration Estimate') +
xlab('Stringency Index'))
## Warning: Removed 36 rows containing missing values (geom_point).
#### Mobility Transit
em <- summary(emmeans(m, ~ Mobility_Transit, at = list(Mobility_Transit = seq(-100,25,by = 10))), type = 'response')
(fig1c <- ggplot(aes(x = Mobility_Transit,y = rEstimatedDuration, col = Session), data = tostat) +
geom_jitter(alpha = .3, size = 1) +
geom_line(aes(y=response), col = 'black', data = em) +
geom_ribbon(aes(y=response, ymin = lower.CL, ymax = upper.CL), col = NA, alpha = .2, data = em) +
scale_y_log10( limits = c(.3,3)) +
scale_fill_manual(values = Palette_Session) +
scale_color_manual(values = Palette_Session) +
ylab('Relative Duration Estimate') +
xlab('Mobility Index'))
## Warning: Removed 36 rows containing missing values (geom_point).
Other mobility measures have little effect on rEstimatedDuration
m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Mobility_Retail + Mobility_Parks + Mobility_WorkPlaces + Mobility_Residential + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
##
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit +
## Mobility_Retail + Mobility_Parks + Mobility_WorkPlaces +
## Mobility_Residential + Subjective_Confinement + Age + poly(Hour_Of_Day,
## 3), data = tostat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.20132 -0.10766 0.00957 0.12317 1.72235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.222e-02 3.418e-02 1.528 0.1267
## Stringency_Index -9.978e-04 4.844e-04 -2.060 0.0395 *
## Mobility_Transit -1.331e-03 7.323e-04 -1.817 0.0694 .
## Mobility_Retail 1.069e-03 7.004e-04 1.527 0.1270
## Mobility_Parks -3.610e-05 2.142e-04 -0.169 0.8661
## Mobility_WorkPlaces 1.122e-03 8.110e-04 1.383 0.1667
## Mobility_Residential 3.119e-03 1.845e-03 1.690 0.0911 .
## Subjective_Confinement -1.217e-03 1.306e-03 -0.931 0.3518
## Age -8.432e-05 2.776e-04 -0.304 0.7613
## poly(Hour_Of_Day, 3)1 -1.612e-01 2.084e-01 -0.774 0.4392
## poly(Hour_Of_Day, 3)2 -2.783e-01 2.055e-01 -1.355 0.1757
## poly(Hour_Of_Day, 3)3 3.045e-01 2.076e-01 1.467 0.1425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1967 on 2135 degrees of freedom
## (267 observations deleted due to missingness)
## Multiple R-squared: 0.01226, Adjusted R-squared: 0.007167
## F-statistic: 2.408 on 11 and 2135 DF, p-value: 0.005686
anova(m)
effectsize::eta_squared(m)
plot(m)
library(patchwork)
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:cowplot':
##
## align_plots
layout <- "
A
B
B"
fig1a / (fig1b / fig1c + plot_layout(guides = 'collect')) + plot_layout(design = layout) + plot_annotation(tag_levels = 'a') & theme(plot.tag = element_text(face = 'bold'))
## Warning: Removed 36 rows containing missing values (geom_point).
## Removed 36 rows containing missing values (geom_point).
imagefile = 'PaperFigures/Fig1.pdf'
ggsave(
basename(imagefile),
plot = last_plot(),
device = 'pdf',
width = 89,
height = 150,
path = dirname(imagefile),
units = 'mm',
bg = "white"
)
## Warning: Removed 36 rows containing missing values (geom_point).
## Removed 36 rows containing missing values (geom_point).